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2D Fokker-Planck models of Rotating clusters 
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ABSTRACT 

Globular clusters rotate significantly, and with the increasing amount of detailed mor- 
phological and kinematical data obtained in recent years on galactic globular clusters 
many interesting features show up. We show how our theoretical evolutionary mod- 
els of rotating clusters can be used to obtain fits, which at least properly model the 
overall rotation and its implied kinematics in full 2D detail (dispersions, rotation ve- 
locities). Our simplified equal mass axisymmetric rotating model provides detailed 
two-dimensional kinematical and morphological data for star clusters. The degree of 
rotation is not dominant in energy, but also non-negligible for the phase space dis- 
tribution function, shape and kinematics of clusters. Therefore the models are well 
applicable for galactic globular clusters. Since previously published papers on that 
matter by us made it difficult to do detailed comparisons with observations we pro- 
vide a much more comprehensive and easy-to-use set of data here, which uses as entries 
dynamical age and flattening of observed cluster and then offers a limited range of 
applicable models in full detail. The method, data structure and some exemplary com- 
parison with observations are presented. Future work will improve modelling and data 
base to take a central black hole, a mass spectrum and stellar evolution into account. 
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1 INTRODUCTION 



Many globular clusters do show so me amount of rotation, 
even the old ones of our own g alaxy jMevlan fc MavorllQSfit 
iLupton. Gunn fc Griffinlll987l and some ten others, K. Geb- 
hardt, personal communication). The more accurate and 
star-by-star observations become available, the more impor- 
tant it is to properly model the effects of rotation together 
with evolution of globular star clusters, since typically the 
amount of rotational energy in clusters is not dominant, 
but also not negligible. It is long known that indeed flat- 
tenings of galactic globular clusters correlates with their ro- 
tation, suggesting tha t rota tion still is important for the 
shape Iwiiitc & Shawl j^S^)- Proper Motions of stars in uj 
Cen have been measured to get a real 3D profile of stellar 
motions including rotation, and in fantastic improvements 
of observations using the Hubble Space Telescope it became 
even possible to measure rotation of 47 Tuc in the plane 
of the sky ijAndcrson & King 2003). While without rota- 
tion the wealth of observational data such as luminosity 
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functions and derived mass functions, color-magnitude dia- 
grams, and population and kinematical analysis, obtained by 
e.g. the Hubble Space Tele scope, for ex tragalact ic as well as 
Milky Way clusters (cf. e .g.lPiotto et alj^20p2: Barm by et alJ 
[ 20021: iPulone et alJl2003l: iRich et al.ll2005l: IZhao et alJl2005l : 



Becc^^^an]200€ , to mention only the few most recently 
appeared papers), is balanced by a decent amount of de- 
tailed modelling (see below), theorists and observers alike 
seem to be appallingly puzzled by rotation in clusters, and 
are disappointed that (multi- mass) King model fitting does 
not work very well (sec e.g. "Piotto e t_ai] llQOfll . and for a 
specific example McLaughlin et al. 20q3, therein). 

One n otable exception being the ea rly three-integral 
models by iLupton. Gunn fc Griffinl (^^3) nobody seemed 
to care about rotation in globular cluster modelling since 
then for a long time. There was a thesis by Jeremy Good- 
man l)l983fl . whose part on rotation of clusters remained 
unpublished, but who was the pioneer for the main idea 
to use a two-dimensional approximation using the distri- 
bution function as function of energy E and z-component 
of angular momentum J^, and neglect the possible depen- 
dencies on a third integral J3 for the time being. Unfortu- 
nately the solution of the 2D orbit averaged Fokker-Planck 
(henceforth FP) equation, including the self-consistent nu- 
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merical determination of diffusion coefficients, is an ex- 
tremely complex numerical problem, involving two levels of 
numerical integrations in a second order integro-differential 
equations, in other words we need 4D nested discretiza- 
tion in E, Jz, vj, and z where the latter are the spa- 
tial coordinates. It turned out that the problem was too 
challenging in the mid-eighties for the then available com- 
puters. With the help of Goodman, wh o kindly made his 
PhD thesis av ailable at an early stage, lEinsel fc Spurzeml 
Paper J) could re-write from scratch a new FP code 
and publish a sequence of rotating cluster models (pre- 
collapse, equal masses). These studies have been r ecently im- 
proved in a comprehensive effo rt on the FP side (iKim et alJ 
I2OO2 Paper 15 12004 Paper 111 , and also promising results 
exist which show that the neglect of the third integral 
does not harm too much (if ffattening is not extreme), 
since there is fair ag reement wit h direct A^-body model s 
(see some examples ii ilPaper iJandlBoilv fc Spu rzcrn '200(f). 
In a new approach iArdi. Spurzem fc Mineshiee ,.2005'1 we 
are working to back up the amount of good Af-body data 
for rotating clusters by a large number of small A'' aver- 
aged models and a fe w large models, Ei s was done fo r the 
non-rotating case bv iGiersz fc Heggi3 il994bl Il996l) and 
iGiersz fc Spurzeml (Il994l) . It is clear that while all work 
known to the authors at this moment concentrates on self- 
gravitating star clusters, the improvement of our knowledge 
and methods in the field of rotating dense stellar systems is 
extremely important for galactic nuclei, too, where a central 
star-accreting black hole comes i nto the game (some sta- 
tionary modelling exis t s, such as iQuinla n fc Shapirolll990l 
iMurphv B. et alJll99lL iFreitag fc Benj|2002^ . Also, claims 
have been made that some globular clusters in the halo of 
our galaxy and aro und M31 contain massive black holes 
iGer ssen et aL|200^, and t hese results have been challenged 
bvfeaumgard tet alJ l|2003f). However, the latest paper could 
only provide a very limited number of case studies, due to 
the enormous computing time needed even on the GRAPE 
computers. Therefore it is very urgent to develop reliable ap- 
proximate models of rotating star clusters with black hole, 
which is subject of current work. 

Dynamical modelling of globular clusters and other col- 
lisional stellar systems (like galactic nuclei, rich open clus- 
ters, and rich galaxy clusters) still poses a considerable chal- 
lenge for both theory and computational requirements (in 
hardware and software). On the theoretical side the validity 
of certain assumptions used in statistical modelling based 
on the FP and other approximations is still poorly known. 
Stochastic noise in a discrete A'^-body system and the impos- 
sibility to directly model realistic particle numbers with the 
presently available hardware, are a considerable challenge 
for the computational side. 

A large amount of individual pairwise forces between 
particles needs to be explicitly calculated to properly fol- 
low relaxation effects based on the cumulative effects of 
small angle gravitative encounters. Therefore modelling of 
globular star clusters over their entire lifetime requires 
a star-by-star simulation approach, which is highly ac- 
curate in following stellar two- and many-body encoun- 
ters. Distant two-body encounters, the backbone of quasi- 
stationary relaxation processes, have to be treated prop- 
erly by the integration method directly, while close encoun- 
ters, in order to avoid truncation errors and downgrad- 



ing of the overall performance, ar e treated in relative and 
specially regularized coordinat es iKustaanheimo fc Stiefe^ 
ll96,4lMikkola fc Aarseth " 1998V their centers of masses be- 
ing used in the main integrator. The factual world stan- 
dard of such codes has been set b y Aarseth and hi s code s 
NBODYn-codes (0 ^ n ^ 7, see lAarsethlll999allhl Eooi) . 
Anot her integrator KiRA aim e d at high-precision heis been 
used (McMillan & Hut 1996'; 'Portcgi es Zwart et alJ Il998l: 
pTakahashi fc Po rtegics Zwart .2000.') . It uses a hierarchical 
binary tree to handle compact subsystems instead of regu- 
larization. 

Unfortunately, the direct simulation of most dense 
stellar systems with star-by-star modelling is not yet 
possible, although recent years have se en a significant 
progress in both har d ware and softw are iMaki no fc 
119981: lAarseth|[l999allbl: ISpurzem|[l99gl: iMakind l2005ll . De- 
spite such progress in hardware and software even the 
largest useful direct A'^-body models for both globular 
cluster and galactic nuclei evolution have still not yet 
reached the realistic particle numbers (A' ~ 5 x 10'' for 
globular clusters, A^ ~ 10^ to 10^ f or galactic nuclei). 
Ho wever, recent work by iBaum gardt. Hut fc Heggie (2002) 
and iBaumgardt fc Makinoril2003l^ has pushed the limits of 
present direct modelling for the first time to some 10^ us- 
ing either NBODY6-I-I- on parallel computers or NBODY4 
on GRAPE-6 special purpose hardware. There is a notable 
exception of a direct one-million body problem with acen- 
tral binary black hole tackled bv iHemsendorf et al.l i2003l) . 
Their code, however, being very innovative for the large A^ 
force calculation, still lacks the fine ingredients of treating 
close encounters between stars and black hole particles from 
the standard A^-body codes. 

Bridging the gap between direct models and the 
most interesting particle numbers in real systems is 
far f rom straightforward, neither by scaling iBaumgardt 
l200lf) . nor by theory. There are two main classes 
of theory: (i) FP models, which are based on the 
direct nu merical solution of the orbit-ave raged FP 
equation fCohn '1979'; 'Murphv R. et alj ll99lY and (if) 
isotropic (Lyndcn-BcU fc Eg glctonTI^Sol: Hcggic 19^^ 
and anisotrop ic gaseous models Jciersz fc Spurzeml 1l993 : 
ISpurzemlll99d) . which can be thought of as a set of moment 
equations of the FP equation. 

On the side of the direct FP models there have been two 
major recent developments. Takahashi (1995, 1996, 199"^ 
has published new FP models for spherically symmetric 
star clusters, based on the numerical solution of the orbit- 
averaged 2D FP equation (solving the FP equation for the 
distribution / — f{E, J^) as a functio n of energy and angu - 
lar momentum, on an [E, J^)-mesh). iDrukier et all (Il999h 
have published results from another 2D FP code based 
on the original Cohn (1979) code. In such 2D FP models 
anisotropy, i.e. the possible difference between radial and 
tangential velocity dispersions in spherical clusters, is taken 
into account. 

Secondly, another 2D FP model has been worked out 
recently for the case of axisymmetric rotating star clusters 
(Paper I; Pap_er_II,) . Here, the distribution function is as- 
sumed to be a function of energy E and the 2-component 
of angular momentum Jz only; a possible dependence of the 
distribution function on a third integral is neglected. As in 
the spherically symmetric case the neglection of an integral 
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of motion is equivalent to the assumption of isotropy, here 
between the velocity dispersions in the meridional plane (r 
and z directions); anisotropy between velocity dispersion in 
the meridional plane and that in the equatorial plane {(f>- 
direction), however, is included. We realize that the evolu- 
tionary models provided by us for rotating globular clusters 
are difficult to use for direct comparisons with observations, 
because they are not easily analytically describable. But 
they are the only ones which fully cope with all observational 
data available nowadays (full 3D velocity data, including ve- 
locity dispersions in zu and 0-direction, rotational velocity, 
density, all as full 2D functions of rzi and z). No other evo- 
lutionary model exists so far which is able to provide this 
information. With t he advent of our new post-collapse and 
multi-mass models JPaoer lit iPaoer IlJ) and the inclusion 
of stellar evolution and binaries (work in progress) we will 
be able to deliver even more interesting resu lts. Alrea dy the 
existing A^'-body study iArdi. Spurzem fc Mineshigd I2065I) 
shows that rotation not only accelerates the coUisional evo- 
lution but also leads to an increasing binary activity in the 
system. 

Here we present a detailed and comprehensive set of 
our model data, covering a fair range of rotation rates and 
initial concentrations of galactic globular clusters, which is 
aimed to enable observers to use them for comparison with 
their data. The basic idea is to start up with a dynamical 
cluster age and a present day observed ffattening, and then 
being able to pick a number of models from our evolutionary 
model data base which allow a prediction (within certain 
ambiguity) of how this clusters started up initially. In this 
paper we present the basic procedure of modelling and how 
we create the data, and give a few exemplary data to show 
the capabilities of our data base which is to be found in the 
web (http : //www. ari .uni-heidelberg.de/clusterdata |. 

Sect. 2 and 3 summarize bricffy the numerical and phys- 
ical procedure used for our evolving rotating cluster mod- 
els, and the choice of initial conditions. Sect. 4 describes the 
comprehensive set of data we provide for every of our models 
using one example set, and Sect. 5 shows how observational 
data can be used to obtain the proper entry points in our 
data for a more detailed theory-observation comparison. 



2 THEORETICAL MODEL 

The stellar system is assumed to be axisymmetric (cylindri- 
cal coordinates {w, z, ip) were used) and in dynamical equi- 
librium. Therefore, the distribution function f{vj, z, ip) is 
almost constant on the dynamical time scale. The only clas- 
sical isolating integrals of a general axisymmetric potential 
(j) are the energy per unit mass: 

£=i„2+0(^,^) (1) 

and the component of angular momentum along the z-axis 
per unit mass: 

Jz = zuv^ (2) 
{E, (j> are negative for bound particles) 

The Boltzmann equation written in terms of E and Jz has 
to be evaluated in the code: 

21] 

dt ) coll 



f{E,Jz,l3) is the distribution function in axial symmetry. 
Here I3 represents a third integral, which is neglected in the 
model (see Paper I for a discussion of possible errors). 

Small angle scattering {Av/v ^ 1) is a condition for the 
expansion of the FP approximation. Diffusion coefficients of 
3. and higher order are neglected. Thus, the collision term 
on the right side of Eq.3 is given by: 

V at y coll V 



-^iiAE)fV)~^^iiAJ.)fV) 



dEdJ, 



i{AEAJ.}fV) 



+ 



1 



2d J, 



r((Ajf)/F) 



(4) 



where V is the volume element given by 2n/vj. 

Due to the relation tr » tdyn (for globular clusters 
is tr ~ lO'"* yr ~ lO^tdyn), the orbit will spread only on a 
relaxation timescale (because of encounters and assuming 
a conserved third integral). Therefore, an orbit-average of 
the FP equation is taken over an area in the meridional 
plane that intersects the hypersurface in phase space with 
the same E and Jz. 

Furthermore, the physical variables E and Jz are trans- 
formed to the suitable dimensionless quantities X and Y. 



The choice of X{E) = ln(j 



where Eq is the en- 



21 d^dl 
dt dt dE 



(3) 



y2c/>a-Ea-E y 

ergy of a circular orbit at the core radius of the cluster, 
guarantees a good resolution in energy space both for small 
energies and for the energy belonging to the central particles. 
y(Jz,-E) = -j^fe- is a dimensionless ^-component of angular 
momentum (0 < \Y\ < 1). A rectangular mesh is generated 
on which the FP equation is discretized. The models pre- 
sented use a number grid size in {X,Y) of (100 x 61). The 
relative error in energy conservation using this configuration 
is less than 0.7 per cent. 

The orbit-averaged flux coefficients are derived from 
the local diffusion coefficients, which are found for the ax- 
ial symmetric geometry via the prescriptions of iRosenbluthI 
(1957) involving covariant derivatives instead of the proce- 
dure employed by Goodman (1983) using a non-covariant 
form. The derivation of the diffusion coefficients makes it 
necessary to specify a background distribution function by 
which test stars are scattered. Using the foreground distribu- 
tion for the background self-consistently would require com- 
putational time proportional to N^^ x Nz x A| x N]^ , where 
the Ni are grid sizes in coordinate and velocity space, which 
would reduce the grids employed to inaccurate descriptions 
of the current problem; thus, as in all previous applications 
concerning 2-dimensional Fokker-Planck methods for stellar 
systems, we set up an appro priate form of the background 
distribution. Herein we follow lOoodmanl (|l983h . giving a ro- 
tating Maxwellian velocity distribution to the background, 
thereby reducing the computational efforts necessary. 

The following steps are taken in the code^: 

• Construction of initial rotating King-models (initial po- 
tential 0(ci7, z) and density n{zu, z) pairs are computed) 

^ Thi s meth od, first developed for spherical systems in ID by 
ICohnI Jl97flh . was originally used bv .Goodman. (.198^ in his un- 
published PhD thesis. 
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• The evolution due to stellar collisions is calculated in 
the F P Step using the method described by iHenvev et alJ 
1)195?/ ) ■ while (f>{vj,z) is constant 

• Recalculation of / due to slow changes in the potential 
(VLASOV-Step). 



3 INITIAL CONDITIONS 

The initial models are rotating King Models of the form: 
/rk(£, Jz) = const X (exp(-/3£) - 1) x exp(-/3J7o J^) (5) 
Initial parameters are: 

oJo ~ •\/9/(47rGnc) • fio dimensionless angular velocity 
Wo — —I3{(f> — 4>t) dimensionless central potential, 
where (3 = CTc is central 1-dimensional velocity dis- 

persion and Uc the central density 

Time units are defined as follows: 
Dynamical time to ICohn 197a) : 



to = 



47rr iV 



(6) 



Where, Mi{= 1) is the initial total cluster mass; rc-{= 1), 
the initial core radius; A'^ is the total number of particles and 
r = 47rG^(Mi/iV)2lnA (In A is the Coulomb logarithm). 
The gravitational constant G is set to 1. 

The initial half mass relaxation time is defined as 
JSpitzer fc Hartll'l97ll) 



= 0.138A 



''hi" 



GMi In A 



(7) 



where rhj is the initial half mass radius. Following the con- 
dition to = 1, evolution of trh can be calculated as: 



U = 2.2087rVA'/r-h3 



(8) 



where M and rh are the current total mass and half-mass 
radius. The cluster radius is given by 



9crc 



AnGric 



(9) 



The dyna mical ellipt i city e dvn of the cluster is calcu- 
lated following iGoodmanI lll983h . as 



STrot + STo-^ — (1 — 2s'^) arccos s — 3sVl — 



(10) 



where s = b/ a = 1 — edyn, Trot is the rotational energy, T^^ 
is the energy contained in the azimuthal component of the 
velocity dispersion and is the energy of all components 
of the velocity dispersion. 

The cluster orbits its parent galaxy in a spherical po- 
tential, at a constant distance (circular orbit and constant 
mean density within the tidal radius). Stars are removed 
instantaneously beyond rt using an energy criteria, i.e. if 
they have an energy larger than Etirt), evaluated at every 
Poisson step. 



4 NUMERICAL RESULTS 

We describe in the Appendix a set of numerical simulations 
classified by ellipticity edyn and age of the system t/trh- The 
database contains following evolutionary parameters: 



• 2D parameters: 

- Distribution function f{E, Jz) 

- Density n{zu, z) 

- 1-dim velocity dispersions a^{zAj,z) (azimuthal), and 
o'i{=az) (in meridional plane) 

- Rotational velocity Vrot{^,z) and angular velocity 

z) 

- Potential (jyivj, z) and anisotropy A{zo, z) 

• Global structure parameters like concentration 
c=log(rt/rc), (rc/rh), escape energy and total mass 

The following is a sample of the database. 

In Fig. we show the dependence of the distribution 
function on the angular momentum for the initial model 
Wo = G, uoo = 0.9 (moderate concentration, high initial 
rotation) at initial time. f{E,Jz) for different Jz and the 
same energy E is represented by lines. Lines of higher en- 
ergy (in the core) are narrower (due to the lower values of 
maximal angular momentum). The slope of each curve rep- 
resents the exponent in the King distribution (Eq. |^ . For 
comparison, we show f{E, Jz) for the same model, at a later 
time (Fig.|5J. We see that the slopes of the curves (rotation) 
change depending on the energy. They are not King-models 
any more, and the initial relation between / and Qo is not 
conserved. Changes in anisotropy (see definition in the Ap- 
pendix) can be also observed in Fig. |5| where orbits lower 
energy (bottom curves) are slightly depleted around J^ — 
(more circular orbits, negative anisotropy). At the same time 
one can observe positive anisotropy built by orbits of higher 
energy (top curves), where / is enhanced around Jz — 
(more radial orbits). 

Although the main source of numerical error was de- 
tected in the recalculation of / due to small changes in 
the potential (VLASOV-step), where a second-order inter- 
polation of / is carried on; the relative error in /, accu- 
mulated over all grid cells per time step, is not larger than 
10~^, which makes the numerical integration stable. The 
sharp bends at both ends of the lines of constant energy 
in Figs. and |5| visualize these errors. Moreover, rotat- 
ing and non-rotating models show nearly the same energy 
conservation accumulated errors (~ 0.7 per cent for a grid 
(A^E,A'jJ=(100,61)). 

Fig. 13 shows the evolution of the central density for 
11 models with different rotation ujq and potential Wo pa- 
rameters . The models with smaller parameter ojq reach col- 
lapse in longer times. Stronger rotating models show smaller 
collapse times. For comparison w e in clude here the models 
Wo=6,7,8. See iPaoer it iPaoer III and iPaoer Il3 for further 
discussion on collapse times. Similarly, in Fig. 2] we show 
the evolution of dynamical ellipticity for the sam e models, 
as presented in lEinsel fc SpurzemI il999. Paper J) (their Fig. 
7). As observed, ellipticity decreases steeply for the initially 
most strongly rotating models and the final states of all 
models lack significant fiattening. Note that the curves can 
cross each other, due to less effective angular momentum 
transport beyond the tidal boundary, in the more moder- 
ately rotating models. In Figs. and |1| evolution is given in 
units of initial trh- 



See Appendix for description of data and access to database 



2D Fokker-Planck models of Rotating clusters - or why there is more than King models 5 

Table 1. Globular cluster observational data (see text for description) 







[ctrcmin] 


T^t [s^rcmin] 




6 


t 1, fvrl 
^rn [J ^ J 




NGC 104 {47Tuc) 


0.40 


2.79 


42.86 


2.03 


0.09 


3.01x10'-' 


1.2±0.1xl0^'' (4) 


NGC 2808 


0.26 


0.76 


15.55 


1.77 


0.12 


1.3x10^ 


1.6±2xl0^° (1) 


NGC 5139 (uj Cen) 


2.58 


4.80 


44.83 


1.24 


0.12 


LOxlO^" 


1.6±0.3xl0i° (7) 


NGC 5286 


0.29 


0.69 


8.36 


1.46 


0.12 


1x10^ 


1.6±0.2xl0i° (5) 


NGC 5904 {M5) 


0.42 


2.11 


28.4 


1.83 


0.14 


3.38x10^ 


13.5±lxl0'J (6) 


NGC 6093 (M80) 


0.15 


0.65 


13.28 


1.95 


0.0 


7.2x10* 


1.4 X 10^" (2) 


NGC 6121 (M4) 


0.83 


3.65 


32.49 


1.59 


0.0 


6.6x10* 


13.5±lxl0S (6) 


NGC7078 (M15) 


0.07 


1.06 


21.50 


2.5 


0.05 


2.23x10^ 


13±3xl0'' (3) 



Refe ren ces.- (1) .Alcaino et al. 1990; (2) Brocato ct al, 1998; (3) iCaputo et al...l984t (i) Icrundahl ct al.. 
I2OO2I: ('5')ISamus et alJll995t (61 Isandauist et al.lll996t (71 Kinematical data from lTrager S. et al.lll995L age 
from lThompson et al.ll200ll 
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Figure 1. Distribution function f{E, Jz) against Jz for the model 
(Wo = 6,uJo = 0.9) at t = Otrhi 



Figure 3. Central density Uc against time t/ty-i^. for models with 
Wo = 6 (dashed lines), Wo = 7 (dotted lines) and Wo = 8 (solid 
lines). The rotation parameter a;o=0.0,0.3,0.5,0.7 (for Wo = 6), 
1.^0=0.0,0.3,0.5 (for Wo = 7) and a;o=0.1,0.15,0.2,0.3 (for Wo = 8) 
increases from right to left at the upper time axis. 



-0.0 0.1 



t/l,„i = 5.534 




Figure 2. Distribution function f{E, Jz) against Jz for the model 
{Wo = 6, LJo = 0.9) at t = 5.53trhi 



5 OBSERVATIONAL DATA COMPARISON 

Observations sho w that flatten i nR is a common feature of 
globular clusters. iGever et al.l lll983h observed the radial 
variation of the ellipticity of 20 Galactic globular clusters 
and 4 clusters in the LMC. They found a general trend 
in ellipticity in a sense that in the inner and outer cluster 
regions they appear more circular and the ellipticity reaches 
a maximum at ~ 1/3 of the observable cluster extension. 
Evhitc & Shawl (1983) have measured the projected axial 
ratio (b/a) of 100 globular clusters in the Milky Way, 
where a and b denoted the semi-major and semi-minor 
radii, respectively. They obtained that the mean axial ratio 
< a/b >= 0.93 ± 0.01 corresponding to < e >= 0.07 ± 0.01, 
where e = 1 — b/a. They argued that the flattened shape of 
the clusters can be caused by either anisotropy in velocity 
dispersion or rota t ion. 

iMevlan fc MavoJ lll98d) and iMerritt et aP (|l993) studied 
globular clusters in our galaxy to reveal without ambiguity 
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Figure 4. Dynamical ellipticity edyn against time t/trh; for mod- 
els witli Wo = 6 (dashed lines), Wo = 7 (dotted lines) and Wo = 
8 (solid lines). The rotation parameter a;o=0.0,0.3,0.5,0.7 (for 
Wo = 6), a;o=0.0,0.3,0.5 (for Wo = 7) and a.;o=0.1, 0.15,0.2,0.3 
(for Wo = 8) increases from bottom to top at the 

^dyn axis. 



Figure 5. Time evolution (in units of current t^h) of dynamical el- 
lipticity. The rotation parameter aJo=0.0,0.3,0.5,0.7 (for Wo = 6), 
1.^0=0.0,0.3,0.5 (for Wo = 7) and aJo=0.1, 0.15,0.2,0.3 (for Wo = 8) 
increases from bottom to top at the e^y^^ axis. 



the global rotation in clusters such as uCen and 47Tuc. 
HST observations of globular clusters provide data on 
the velocity field o f glob u lar clusters to high acc uracy 
JVan Leeuwen et al.l I2OO0I: lAnderson fc Kind l2003l) . We 
look forward to have direct integration models to compare 
with these data. 

Table gives observational data for 6 globular clusters 
with ellipticities larger than zero and 2 with e = 0. The data 
was take n from the H arris Catalog of Milky Way Globular 
Clusters lHarrislll99el. February 2003 revision). In the table 
are: 

Cluster identification number 

current core radius in arc-minutes; 
current half mass radius in arc-minutes; 
current tidal radius in arc-minutes; 
logjQ(rt/rc), concentration; 



rt, 

C 

e, ellipticity (projected axial ratio); 
trh, half mass relaxation time in years; 
cluster age in years (taken from the literature, as 



Column 1: 
Column 2; 
Column 3: 
Column 4: 
Column 5: 
Column 6: 
Column 7: 
Column . 
remarked) . 

Note that the listed values of c and rc of M15 should 
not be used to calculate a value of tidal radius rt for this 
core-collapsed cluster as explained in the catalog. 

A first approximation is achieved by using the age of the 
system in units of current half-mass relaxation time (t/trh) 
and the cluster ellipticity. Fig. |3 shows the dependence of 
ellipticity on the cluster age, where the position of the iden- 
tified models for six clusters of Table are included. For 
the model identification, the errors in the cluster age de- 
termination (taken from the literature. Table define age 
intervals, which match with the co mputed el l ipticity cu rves. 
As commented in previous papers ijPaper l lPaoer Iljl . the 
observed decrease in ellipticity is not only due to angular 
momentum loss (evaporation) but also due to the expan- 
sion of mass shells retaining their angular momentum and 
decreasing their angular velocity (inversely to the actual ra- 
dius of the shell) . A coincidence in the behavior of ellipticity 
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Figure 6. Evolution of concentration. The rotation parameter 
0.^0=0.0,0.3,0.5,0.7 (for Wo = 6), aJo=0.0,0.3,0.5 (for Wo = 7) and 
0^0=0.1,0.15,0.2,0.3 (for Wo = 8) increases from top to bottom at 
the concentration axis. 



and rotational veloc ity profiles gives an evid ence of rotation 
in globular clusters ijMevlan fc Mavo3ll98fih . 

As far as concentration and ellipticity are independently 
computed, we use the age and concentration c = log(rt/rc) 
of the cluster, in order to make the determination of the 
model more accurate. Fig.|S|shows the concentration against 
dynamical age for the same models. Using Figs. and |S1 
one can determine the initial parameters (Wo, wo) of the 
model to use. 

Table |5| shows the cluster parameters for GCs of Ta- 
ble Q with e > 0, used to determine the initial conditions 
of the simulation (Wo,wo). First column shows the name of 
the cluster, second column the observed ellipticity, third col- 
umn the current concentration, fourth column the age of the 
cluster in units of half-mass relaxation time and columns 5 
and 6 the initial parameters Wo (King parameter) and uiq 
(initial rotation). 

For the model determination the coincidence on both, 
ellipticity and concentration, to observations in Figs. |^ and 
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Table 2. Globular cluster initial models (see text for description) 
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c 


t/trh 


Wo (initial) 


(initial) 


47Tuc 


0.09 


2.03 


4.0±0.3 


8 


0.15 


N2808 


0.12 


1.77 


12±1.5 


8 


0.3 


ujCen 


0.17 


1.24 


1.6±0.3 


6 


0.5 


N5286 


0.12 


1.46 


16±0.2 


7 


0.5 


M5 


0.14 


1.83 


3.99±0.3 


8 


0.2 


M15 


0.05 


2.5 


5.8±1.3 


8 


0.1 



|S|was located. As these figures could be filled using all pos- 
sible pairs {Wo, uq), we show representative initial parame- 
ters and rule out the ones which do not match both initial 
conditions (e.g. higher uiq, or lower Wq). Thus, determined 
models are not unique but representative. Nevertheless, as 
the models are still idealized, implementation of more com- 
plexity in cluster structure (different mass species, stellar 
evolution) and in interaction with the parent galaxy (tidal 
shocks, dynamical friction) could affect the results. 

In Figs. 171 IHl l^and llll we present contour plots of the 
rotational velocity v^ot = ti7 f2 in the meridional plane {zu, z) 
for four models of Table |21 Wrot was transformed to physical 
units (km/s) using an estimation of the total cluste r mass 
according to a mean mass to light ratio of M/L = 2 jHarrij 
Ll996, February 2003 revision) and using Vc as a length scale. 
The distances {zu and z) are rescaled to arc minutes. Single 
mass pre-coUapsed models were used. 

Fig. |7| shows a contour map of rotational velocity 
for an initial model {Wo = 8.0, ojo ~ 0.15) at f/frh = 
4, identified as a 47 Tuc-model. The maxima of v^ot ~ 
7.6 km/s is located at ~ 7.7'. Rotation HST measurements 
of 47 Tuc have recently confirmed earlier studies of this 
globu l ar cluster l[^]^ersor^^ Kinell2003t iMevlan fc MavoJ 
ll986l) . lAnderson fc Kin^ T200^) found a rotational velocity 
of 5.7 km/s at a radius of 7.5 pc. This is quite comparable 
to the maximum line-of -sight rotation of ^ 6 km /s at 5 pc 
(~ 11 - 12 rc) found bvlMevlan fc Mavoil jlQSd) . 

iGerssen et al.l ll2002l) reported the kinematical study 
of the central part of the globular cluster M15, including 
Vrot- M15 is kno wn as a globular clus ter which contains 
a collapsed core. iGebhardt et al.l lll994l) showed that M15 
has a net projected rotation amplitude of v^ot = 2 km/s 
at radii compa r able to the half-light radius (about 1'). 
IGebhardt et al.l i2000h has revealed that the rotation am- 
phtude is larger {v^ot = 10.4 ± 2.7 km/s) at r ^ 3.4", 
implying that v^ot /cr ~ 1 in this region. Fig. |H| shows the 
rotational velocity distribution in the meridional plane of 
a M15-model (W^o = 8.0, ujq = 0.1) at t/t,i, = 5.8. Al- 
though our single-mass (and multi-mass) models can not 
explain the rapidly rotating central region, this phenomena 
was already observed in our current simulations of rotating 
collapsed BH-models, which will be topic of a future publica- 
tion. Note that the core radius of this cluster (0.07') requires 
a better resolution to be able to confirm observational data 
(see below). 

iGerssen et al] J200l l2003l) . used HST to obtain new 
spectra of stars in the central cluster region of M15. They 



47TUC 




Figure 7. Contour map of rotational velocity for a 47Tuc model. 
Wo = 8, ujo = 0.15 and t/tj-h = 4.0. Distances are in arc minutes 
and contours are labeled in km 
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Figure 8. Same as Fig.[7|for a model of M15. Wq ■ 
and i/trh = 5.8 



UJO 



0.1 



found theoretically, that the existence of an intermediate- 
mass black hole of mass (1.7 ± 2.7) x 10"^ A/ can explain 
the observed velocity dispersion. However, IgaumEardt et ^ 
|2003) have shown that the core-collapse profile of a star 
cluster with an unseen concentration of neutron stars and 
heavy mass white dwarfs can explain the observed central 
rise of the mass-to-light ratio. Though, they did not reject 
the possibility of the presence of IMBHs. 

Fig. |n| shows the contour map of rotational velocity 
for an initial model {Wo = 6.0, tJo = 0.5) at i/trh = 1-6 
identified as a wCen-model. The maxima of Vrot (6 km/s) 
is found at ~ 9' ~ 3rc. For comparison. Figure [TI71 shows 
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Figure 9. Same as Fig.lTlfor a model of tjCen. Wo = 6, ojq = 0.5 

and t/trh = 1-6 



Figure 11. Same as Fig. |2|for a model of M5. Wo = 8, ujq 

and t/trh = 4.0 



: 0.2 




Figure 10. Contour plot of the mean azimuthal velocity in the 
meridional plane of ijCen. Distances are in arc minutes and con- 
tours are labelled in km/s (Merritt et al. 1997). 



the observed two dimensional rotation al stru cture of tjCen 

iDliii3)- 



in the meridional plane l|Merritt et al.Ml997l) . Using a non- 
parametric technique on 469 radial velocity data they ob- 
tained a two dimensional structure of rotational velocity in 
wCen and found a peak rotational velocity of 7.9 kms~^ 
at ~ 11 PC correspond i nK to ~ 3 times the core radius 
fr^ 't. lVan Leeuwen et alJ 1)20 00^ found agreements of proper- 
motion studies with the rotation in this cluster. 

Fig. 111! shows the rotation in the meridional plane of a 
M5-model with a maximum of Urot = 5.6 km/sec at about 
5' ~ 12rc. It has one of the highest ellipticities in the Harris 
catalogue (e = 0.14) and the highest in our sample. Cur- 



rently there is not observational data of rotation of this 
globular cluster to compare with. 

The simulated profile of frot/c" (rotational velocity over 
one dimensional velocity dispersion) of four clusters of Ta- 
ble|21is plotted in Fig. 1121 They show a maximum of rotation 
at around the half-mass radius and a rapid decline after- 
wards, as well as a rigid body rotation in the center. Note 
that the max;ima of the M15-model is located closer to the 
center than the other clusters, as M15 has the largest dy- 
namical age (see Table|5J and dynami cal effects as gravogyr o 
instabilities dominate the evolution jHachisul 11971 1198211 . 
The results show that our data i s consistent with observa - 
tions, as the newly presented in IVan de Ven et alJ i200d) . 
who show the importance of rotation (vrot/o" distribution 
in the meridional plane) in ojCen at radii of about 5 to 15 
arcmin, however dropping more inwards and further our- 
wards. Nevertheless, more physical constraints (multimass 
modes and stellar evolution) are needed in our simulations 
in order to follow properly the interaction of rotation and re- 
laxation. In the same way, we can obtain data of any model 
described in Section 4 and in the Appendix. 



6 CONCLUSIONS AND OUTLOOK 

Most globular clusters are flattened by rotation, and show 
interesting features in their rotation curves. While rota- 
tion is not dominant it is also not negligible, the total 
amount of rotational (ordered) kinetic energy (as compared 
to unordered kinetic energy) can be as high as a few 10 
% typically. Young globular clusters are usually more flat- 
tened than older ones. This was observed in the Mag- 
ellanic clouds (Kontizas et a l., ,199Q: Hill fc Zari tskv 20o|) 
and in M31 l|Barrnl^^erST|20(d) ; and theoretically stud - 
ied jFrenk Fallll 982HOnndwinlll 9971: iRekki et AM)^. 
But see also iMackev fc van den Ber3 J2005l) for the MW. 
While observations, especially using the utmost accuracy for 
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Figure 12. Simulated v/cr plotted against cluster radius, given 
in units of current half mass radius. 



single star measurement s in globular clusters with the Hub- 
ble S pace Telescope (e.g. lAnderson fc Kinel2003l : lRich et alJ 
|2005^ provide a marvelous detail of morphological and some- 
times surprising kinematical information, people seem some- 
times still surprised that even the most elementary fea- 
ture of flattening and rotation (and its dynamical impli- 
cation for rotational and dispersion velocities) are so dif- 
ficult to model and will not be fitted well by multi-mass 
King models (Longaretti & Lagoutc 1997; Paper 11^. We 
have provided detailed evolutionary models here which make 
it possible to compare t hem easily with observations, us- 
ing data first obtained bv lEinsel fc Sourzeml (Il999. Paper J) 
for the equal mass pre-coUapse phase. The basic idea of 
our data collection, which we have discussed here in a 
few examples, and whicfi is complete at the web URL 
jhttp : //www. ari .uni-heidelberg . de/clusterdata is that 
one can go in it using simple global observational data (dy- 
namical age, flattening and/or concentration of a cluster), 
and then pick from a small number of theoretical evolved 
models to compare with. 

2-dimensional axisymmetric FP models were used in or- 
der to follow the dynamical evolution of rotating stellar sys- 
tems by using initial parameters (King dimensionless central 
potential parameter Wo, dimensionless measure of rotational 
energy uq)- Our data sample is the only one available from 
theoretical evolutionary modelling of rotating globular clus- 
ters and is aimed to allow direct comparison w i th obs erva- 
tional data, as t he pr esented by iMerritt et alj (^93) and 
iGebhardt etall tOOd) . by using our formatted data files 
and IDL plotting routines provided. Rotation conduces the 
system to a phase of strong contraction (during a mass and 
angular momentum loss process) and a ccelerates the col- 
lapse time (gravo-gyro phase; e.g. lHachisulll979lll98l . The 
observed rotational structure in the meridional plane and 
Vrotlo curves (despite of the central region of clusters like 
M15, discussed in Sect|SJ can be reproduced by our models. 
Nevertheless our models are still idealized and simplified. 
Physical scalings, such as time scales, may change if more 
effects, such of binaries, are included. 

The provided database of models with theoretical clus- 
ter data will make it possible for observers to select model 



data to compare with Our database will in the future 
be extended to more realism including a stellar mass spec- 
trum and stellar evolution of singles and binaries, which all is 
subject of ongoing work CArdi. Spurzcm fc Mjpe shige 200^. 
Also a study of rotating clusters with a central massive star- 
accreting black hole is under way (Fiestas et al. 2006, in 
prep.), which will provide important data for galactic nu- 
clei and to check claims that some globular clusters contain 
intermediate mass black holes. 
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APPENDIX A: DATA DESCRIPTION 
Grid of models 

In http://www.ari.uni-heidelberg.de/clusterdata/ we 
present exemplarily a grid of runs of our 2D FP models 
Wo = 3,6,9 and loo = 0.0 (no-rotating), 0.3,0.6,0.9 
(rotating). The data is classified by ellipticity and age of 
the system. Following, a briefly description of the data. 

The initial classification of the models is showed in Ta- 
ble 1 (on the web). {Wo,ujo) are the respectively King- and 
rotational initial parameters. Second table shows the grid of 
database. Each point of the grid represents an ellipticity and 
time interval. The numbers in the cells give the total number 
of available datasets with the corresponding age-ellipticity 
relation. 

By selecting one grid point, one gets the datasets in- 
cluding following information (see Table lAH : 

(i) Structural parameters: 

concentration c = log(rt/rc), (r'c/rh), escape energy, total 
mass and energy at any time of evolution. 

(ii) density n{m, z) 

(iii) Angular- {Q,{jii,z)) and rotational velocity 

Wrot(l^, z) — Q,{'ClI, z) ■ 

(iv) 1-dim total velocity dispersion at{zu,z), azimuthal 
vel.disp a,p{zu,z) and vel.disp. in meridional plane ar(cc7,z) 

(v) Potential (j}{vj,z) and anisotropy A{zu,z). 

(vi) Normalized (unnormalized) energy X{E) and angu- 
lar momentum Y{Jz), and evaluated distribution function 
f{E,J.). 

Further definitions: 

• escape energy xq = '^^''^^^''^'^ (^(j)^; central potential). 

• core radius rc(Eq.|nil 

• 1-dim square velocity dispersion: at{uj,z)^ = 
a^[TXj, z)'^ -\- 2(Jr(t47, z)^, (due to ar — <j^) 

. Anisotropy = 1 - ^^g^ 

• Normahzed energy X{E) = ln ^'^'=~|~-^° 

• Normalized angular momentum Y{Jz,E) = j^'I^^^^ 

Evolution in time 

Here we present filtered data classified by initial model 
parameters {Wq,u!q) (Table 3 on the web). They include 
the time evolution of following parameters: 
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Table Al. 













File name 


brief description 


GP.[Wo 


X 


10]. 




X 


100].[timcstcp] 


global parameters 


RHO.[Tyo 


X 


10]. 


[<^o 


X 


100].[timcstcp] 


density 


ROTV.[Vyo 


X 


10]. 




X 


100].[timestep] 


rotational- and angular velocity 


VDisp.[vyo 


X 


10]. 


[<^0 


X 


100].[timestcp] 


ID velocity dispersion and components 


PHI.AN.[Wo 


X 


10]. 




X 


100].[timcstcp] 


potential and anisotropy 


DF.[Wo 


X 


10]. 


[<^0 


X 


100].[timcstcp] 


distribution function 



• time in units of half mass relaxation time (t/trhi) 

• time step number 

• time in units of central relaxation time (t/tc) 

• core radius rc (units of initial core radius=l) 

• central density Uc 

• central rotational velocity Vc 

• central 1-dim velocity dispersion ctc 

• ellipticity edyn (Eq. 1101 

• escape energy 

• collapse rate ^ = ^ ^ 

• total mass (units of initial) 

• total ang. momentum (z-component) 

• total potential energy 

• total kinetic energy 

Plots as those presented in this Paper are also shown in- 
cluding IDL-routines used to generate them. This paper has 



been typeset from a TJriX/ M^jX file prepared by the author. 



